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Abstract 

We study the Schrodinger equation which comes from the paraxial approximation of the 
Helmholtz equation in the case where the direction of propagation is tilted with respect to 
the boundary of the domain. In a first part, a mathematical analysis is made which leads to 
an analytical formula of the solution in the simple case where the refraction index and the 
absorption coefficients are constant. Afterwards, we propose a numerical method for solving the 
initial problem which uses the previous analytical expression. Numerical results are presented. 
We also sketch an extension to a time dependant model which is relevant for laser plasma 
interaction. 

1 Introduction 

For the simulation of the propagation of a monochromatic laser beam in a medium where the 
local refractive index is nearby a constant, it is classical to use the paraxial approximation of the 
Maxwell equations. This approximation takes into account diffraction and refraction phenomena ; 
it is intensively used for decades in optics and in a lot of models related to laser-plasma interaction 
in Inertial Confinement Fusion experiments (cf |1],[I0], [21], [H] and the bibliography of these 
references). Let us first recall briefly the outlines of this approximation. Denote by 27re the laser 
wave-length, it is in the order of 1 iim and is very small compared to the characteristic length of 
the simulation domain (which is in the order of some mm for the Inertial Confinement plasmas). 
According to laws of optics, the laser electromagnetic field may be modeled by the solution ^p of the 
following Helmholtz equation (which comes from the time envelope of the full Maxwell equations): 

e'^Aip + V + 2ievt'ip = 0, (1) 

where we have denoted: 



ft(x) = z^(x) + i/x(x), 
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so ut is a complex function, its real part v corresponds to a conveniently scaled absorption coefficient 
and its imaginary part to the variation of the refractive index (1 — 2e/i is equal to the square of 
the refractive index n up to a multiplicative constant). 

We assume also that the light propagates according a fixed direction defined by the unit vector 
k. After making the classical WKB expansion: 

V' = ^/exp(— — ), (2) 

equation ([T|) may read as 2ivtu + 2ik.Vii + eA_Ln = e(k.V)^n, where A_l is the Laplace operator 
with respect to the transverse variable: 

A_|_» = V.[(l — k (g) k)V»], 1 being the unit diagonal tensor. 

Assuming that u is slowly varying with respect to the longitudinal variable, we can neglect the right 
hand side of the previous equation. Therefore u satisfies the classical paraxial equation for wave 
propagation: 

zk.Vn + — A_|_n + ivtu = 0, with vt = v + i/i. (3) 

For this kind of model, it is usual to handle a simulation box which is a parallelepiped and 
the laser beam is assumed to enter into the simulation box on a plane boundary denoted by Tq. 
Let us denote n the outward normal vector to the incoming boundary Tq. Classically, the crucial 
assumption is that the laser beam enters into the simulation domain with a very small incidence 
angle, that is to say the vector k is almost equal to — n. Then, in such a framework, ([31) is a classical 
linear Schrodinger equation, the operator k.V plays the part of time derivative and the boundary 
condition on Tq which reads u = u^^ (where u*" is a given function defined on Tq) plays the part of 
the initial condition. On the other hand, artificial absorbing boundary conditions are to be imposed 
on the faces of the simulation domain parallel to the vector k, (see for example [I], [7j, [E]). The 
numerical methods are always implemented on an orthogonal mesh and are based on a splitting 
with respect to the main spatial variable between the diffraction part (|Axn) and refraction part 
{wtu), see [3], ^0] for example. 

We address in this paper a different case where the incidence angle of k with — n is large; these 
simulations are called tilted frame simulations. This kind of simulations is of particular interest 
if one has to deal with the crossing between two beams (in the high energy laser devices, a large 
number of beams are focused on the target, therefore beam crossing may be taken into account, see 
[8] for a survey on related laser propagation problems); an example of such simulations in a very 
simplified case may be found on Figure [131 This tilted frame model has been considered some years 
ago by physicists for dealing with beam crossing problems (see [20|). 

Simulations in a tilted frame are also necessary for dealing with special situations. For instance 
for the propagation of a beam in a domain where the profile of the refractive index n is such that 
n^(x) = nQ(l — e/i(x)) (with no constant smaller than 1) in a first subdomain V and n^(x) = 
AA(x.n*) + (5AA(x) (where J\f G [0,no] depends on a one-dimension variable x.n* and 6J\f is small 
with respect to 1) in a second juxtaposed subdomain P^, one must handle the paraxial equation 
jS]) in subdomain V and the Helmholtz equation ([T|) in subdomain . For the numerical solution 
of ([1]), one has to solve a huge linear system (corresponding to the discretization of the equation 
on a very fine grid) and for handling this huge linear system, it is necessary that the variable x.n* 
corresponds to one of the main direction of T)^ . Therefore the full simulation on (T> U T>^) has to 
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be performed in a box such that the corresponding normal vector n must be parallel to n* (see [6J 
for details for this kind of simulations). 

In the case of a large incidence angle, the crude expansion ^p = f7exp(— m.x/e) leads to dif- 
ficulties and to overcome these difficulties, it has been proposed in |i3j to replace the transverse 
Laplacian by a pseudodifferential operator, but with this approximation, U is not slowly varying 
with respect to the spatial coordinates therefore it is necessary to handle very fine mesh -at least 
10 cells per wave length- to get accurate results. One can also refer to the works in the spirit of [E] 
in the acoustic framework but the application to the optics problems seems to be difficult. 

Here we consider the expansion ^p = uexp(zk.x/e), with u slowly varying with respect to k.x, so 
we have to deal with the tilted frame Laplace operator Aj^ and one has to supplement the equation 
([3|) with a right incoming boundary condition on Tq. For the statement of this boundary condition, 
one assumes that a fixed plane wave = u*" exp(ik.x/e) enters into the domain where u*" is 
a given function of the variable which is orthogonal to k. Now, for the Helmholtz problem, the 
boundary condition is classical and may be written as (en.V + ik.n)('i/' — u'^^e^^'^^^) = 0, then using 
([2|) and an asymptotic expansion with respect to the small parameter e, the corresponding boundary 
condition for equation jS]) may read in a natural way as: 

(en.Vx + 2ik.n)(n - u^") = 0, (4) 

where = V — k(k.V) denotes the gradient orthogonal to k. See [9J for a justification of the 
paraxial approximation in the special case we are dealing with. 

If one sets x = {x,y,z) in 3D and x = (x,y) in 2D, the entrance boundary Tq corresponds in 
this paper to x = 0. In the sequel we consider a 2D problem but most of the ideas of this work may 
be extended to the 3D case. 

Equation ([3|) may be recast as: 

i{kxOxU + kydyu) + + ^'^t'^ ~ 

and up to our knowledge, the numerical solution of this kind of equations is novel; the main difficulty 
is to handle correctly the tilted Laplace operator A±u. For the mathematical analysis of the 
problem, one key result is the following (cf. proposition 2). On the half-space {{x,y) s.t. x > 0}, 
if the coefficient vt is a positive real constant, after taking the Fourier transform with respect to 
the y variable, the problem ©(III) is equivalent to an ordinary differential equation with respect to 
the X variable and it is possible to exhibit an analytical solution. This analytical formula is the 
convenient tool for numerical treatment of the diffraction part of ([3]) in the general case where vt is 
not constant. 

The paper is organized as follows. In Section El after setting classical energy estimates for 
Problem jS]) supplemented by we prove the above mentioned theoretical result. 

Section [3] is devoted to the description of the numerical scheme for solving Problem jSlldll) ; it is 
based on a splitting method with respect to the spatial variable x using fast Fourier transforms on 
a first step (for the diffraction part) and a standard finite difference method on a second step (for 
the advection and refraction part). 

In Section HI we give the numerical results on the initial problem and for a model where the 
coefficient /i in ([3]) is replaced by /(|ti|) corresponding to the autofocusing which occurs in the 
laser-plasma interaction (see [19] for instance). Prom a physical point of view, this term represents 
a variation of the plasma electronic density caused by the ponderomotrice force of the laser. In 
the last section we consider a more general model where the stationary problem ([3]) is replaced by 
a time dependent one which is coupled to a hydrodynamic system for a suitable modeling of the 
plasma behavior. 
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2 Analysis of the Tilted Paraxial Equation 

For reasons which will appear in the sequel, we assume in this section that 

infxi^(x) > 0. (5) 

We first study the problem where the simulation domain is the half-space: 

P = {x=(x,y) s.t. x>0}, ro = {x=(0,?/)}. 

Assuming that /x is a bounded function, we consider the following problem: 

ik • \7u + —A±u — fiu + iuu = on P, (6) 
(ien.V_L -2k.n)(n-n™) = on Tq. (7) 

2.1 Energy Estimate 

Let us first state the following classical estimate. 

Proposition 1 Let (ien.V_L - 2k.n)tt'" G L2(E). If u £ H^(V) is a solution to Problem © 

it is unique. Moreover, we have the following stability estimate, with a constant C independent of 

2v\u\^ + j |k- nllupdy < C j | (ien.V± - 2k.n)'u™pdj/. 



V 



To 



Proof. Let us denote D = n.Vx. Doing the scalar product of Equation jS]) with u and taking 
its imaginary part, we get: 



+ ^(uZ?n - -uZ^u) )dy + / / 2i/|up(ix = 0. 



To V 
According to the boundary condition ([7]) we check that: 



2i 



{uDu - uDu) = -2k • n\u\^ +lm{u{eD + 2ik • n)n*''). 



Then we get: 



2i/|up(ix + j \k.n\\u\'^dy = -Im{J u{eD + 2ik.n)u'''dy) . 
V To To 



(8) 



According to ([8|), if {ieD — 2k • n)u*" = 0, we see that // 2z/|up(ix = 0, so u = 0. Therefore we get 

V 

the uniqueness of the solution of Problem ([6]) JT]) . 

To obtain the stability inequality, we first see that Equation jS]) implies: 



k • n| / \ur < 



To 



j |m|2 j |(eL> + 2fk-n)M*"|2. 

ro \ro 
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Using this estimate, Equation ([8|) leads to: 



2i/|tiPdx + / Ik • nllnp < 



V 



\u\ 



\ To \ To 



|(eD + 2ik-n)ti™|2 < 



Ik • n| 



|(eL> + 2ik-n)u 



in. 1 2 







By the same technique we get also the following estimate: 
|k.n| I {ieD + 2k.n)U|2 

To 



2i/kr + 



2|k.n| 



, ,, , ,2 l,(ieL>-2k.n ii*",2 

k.n ur + - -^^ — — r 

" ' ' 2' 2 k.n ' 



V To 

which says that the absorbing energy plus the the outgoing energy is equal to the incoming energy. 
2.2 Analytical Form of the Solution in the Case vt Constant 

We now assume that = and v is constant for getting an analytical form of the solution to 
Problem (l3|)(|l|. We denote k = (/c^, ky) and g the function defined by: 



2A^a;5 — i^kyl^kxdy hydx)u + 2kxU . 
The problem may read as: 



i{kxdx + kydy)u + - {kxdyy — 2kxkydxy + kydxx)u + ii^u = 0, on V, 



i€ky{kxdy kydx)u + IkxU — IkxQ, 



onrn 



(9) 

(10) 
(11) 



In the sequel, the Fourier variables related to x and y respectively are ^ and ry. The Fourier transform 
in X and y are denoted by J^xi") and J^y{»), moreover J^y{u;x,.) denotes the Fourier transform of 
u{x, .). 

Here and in the sequel, denotes the principal determination of the square root (its real part 
is positive). Denote: 



y y 



' 1.2 



+ 2fz.-^). 



Since > 0, one can define R- without ambiguity and one checks that lZe{R^{iri)) < for all rj. 
Let 5'(1R) be the space of tempered distributions. 

Proposition 2 Assume that g G 5'(1R), then there exists a unique distribution u{x,.) continuous 
from R'*' into S'y(R), solution to Problem / f70j) / f77]) . It is given by: 



Ty{u\x,ri) 
It satisfies also: 



1 + A 1 



(12) 



dx - R-{iri) \Ty{u]x,ri) = 0. 
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Proof. 

The principle is to take the Fourier transform in y of the problem, and afterwards we shall 
consider Fourier transform in x of the equation extended to the whole space. 

Let u be a solution of Problem (fTO])(fTT]) and v the extension of u by zero in the whole space: 
v{x,y) = u{x,y)lx>o- By introducing formally the function v in Equation ([TO]) we get: 



:0- 



ik-Vv + ^A_iV + ii^v = (^{ikx - ^{2kxdy - kydx))u{0,y)^6x=o + ^u{0,y)5'^. 
The term dxu{0,y) is defined by the entrance boundary condition (fTTl) . so we get: 

ik-Vv + I A_Lf + ii^v = ikxg{y)5x=o - ^ {k^dyU^O, 2/)4=o - kyu{0, y)S'^=o) ■ 

Assuming that u £ C(1R+, we are allowed to take the Fourier transform of this expression. 
Let us define P{X, Y) as the polynomial which characterizes the differential operator of the equation, 
that is to say: 

P{dx, dy) = i{kxdx + kydy) + ^{kpl - 2kxkyd% + k^d^J + il^. 



Writing uo{y) = u{0,y), the Fourier transform in y of the equation in v reads: 

■vy 



P{dx,ir])Ty{v;x,rj) = -^l \ -^^yi9\v) - i^r]J='y{uo;r]) ]5x=o + J^y(uo; ??)4=o 



Polynomial P may be factorized as: 



P{d.,iv) = ^ (dx - R+{iv)^ (dx - (iv) ] , (13) 



where we define R±{iv) = i^V " ± ^1 - 2^ + 2ii/^^ . Thus: 

dx - R+{ir]) ]{dx- R-{iri) ]j^y{v;x,T]) = 



'^^^'^ ^y{9;ri) - ^^??-^y(^io;f?)^'^x=o +-^j/N;??)'5l=o- (14) 



ek'^ ny 



We now show that there is a unique acceptable solution for this ordinary differential equation. Let 
us take its Fourier transform in x: 



- R+{iri) R-{iv) ]:Px^y{v;^,ri) = —r^^yidW) - ^(tt^ ~ O^y{uo;v)- 



bmce we can divide each side of this equation by ■^P{i^, irf) : 



— i?+ (irj) — i?_ {irj) ' 

R-{iri)—i—ri 

where a±(r?) = ± R^(i^)_R_'^,^) -^i/(no; ??) ± ^ R^ii^)-R_(^i^) ^y{9\^)- 
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IfO e C\E, one knows that: 

1 _ r ^,(l,>oe^";0 if7^e(0)<O 



I -•^x(lx<oe^";e) if7^e(0)>O. 

Here TZe{R+) = —TZe{R-) > 0. According to the previous remark, since v{x, •) = for x negative, 
one gets a~^{T]) = and 

^,(u;x,77) =a-(r?)e^-(^'')^l,>o, 
so we get J^y{uQ;r]) = — ^^"^-^^yT^- Equahty (fT2]) and the last assertion follow. 



Notice that we can easily calculate, with this formula, the value of the derivative k • Vu. As 
soon as u is regular enough, we can perform an asymptotic expansion in e and u, and find: k • Vu = 
0{e + u). 

From this result, one deduces the following stability result. 

Corollary 1 If g £ H~2(Ji) then the solution u to Problem / f70|) f77] ) is continuous from into 
L^(IR), and it satisfies, for some constant C not depending on the coefficient v: 

lkllLg°(R+,L2(R)) < C\\g\\^_^^^^. 

Since C does not depend on the absorption coefficient i/, one can check that if n*" is smooth 
enough, for x fixed, the function .) converges strongly to a function in l?y when ^ 0. Therefore, 
one may claim that there exists a bounded solution u to Problem ([TO]) (fTTI) . even if = 0. 

Proof. 

Let us integrate with respect to r\ the square modulus of both sides of Equation (fT2]) . Since 

|gR_(ir?)x| _ ^■Re{R-{iri))x ^ ^ g^^^. 

it suffices to show that there exists a constant Ci > 0, not depending on i^, such that: 



l + \v\^<C- 



4 



Vr/ elR. (15) 



2ek k'^ 

So, if we denote X = 1 r^^ry and A'" = 2ev-^^ one first sees that: 

|1 + VrTIiV|2 = 1 + Vx2 + Ar2 + 2(X2 + N^)\cos{^ - ^^E^^^illL) > yTTx^ 

(indeed the cosine is nonnegative) . With a = we have 1 + |ryp = 1 + a^(l — X)^ and it is easy 
to check that 1 + 0^(1 - X)^ < Ci(l + X'^) for Ci = + 1 ; Inequality ^ follows. 

Remark: with the same techniques, one can also find existence and uniqueness of a solution in 
other spaces, for instance, if ^j^'^^^p^i/s £ ^'^(M)^ we have u £ L'^iV). 

Since \Ty{g]ri)\ < C(l + |77p)-^/^|jrj,(n*"; r/)|, that means that if is smooth enough (in H^^^ 
for example), the solution u belongs to L'^{'D). 
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2.3 Remark on the Problem on the Quadrant 

We now consider the same problem (fTO]l(fTT]l but restricted to the quadrant {{x, y) s.t. x > 0, y > 0}. 
To find a good absorbing boundary condition on the boundary {y = 0}, we formally factorize the 
differential operator of Equation (fTOl) as follows: 

P{d., dy) = e^{dy - A^{d,)) {dy - (l6) 

where A^{.) and A^{.) are the roots of P considered as polynomials in dy : 

The definition of the fractional derivative is classical and is based on Fourier transform. The 
quadrant problem that we consider consists of Equations ([TO]) (fTTl) supplemented with the following 
boundary condition 

dyu- A+{d^){u) ={), Vx>0, for y = 0. (17) 

Then, we have the following result, which is detailed in [HE] (for related boundary value problems 
for classical Schrodinger equations, see for example |12|). 

Proposition 3 Assume g G //~2(]R,+ ) and its support is in (0,+oo). Let u be the solution of the 
half-space problem / fTOjJ / fTT]) . There is a unique solution U continuous from 1R+ into L^(IR^) of 
Problem / flOj) / flT]) / [77|) and it satisfies 
i) if ky > 0, then U = uly>o, 

a) if ky < and if the incoming data is given by g{y) = h{y — a) with a > 0, then: 

aij™oo ~ ^1j/>o||l°°(]r+,l2(r+)) = 0. 

3 Numerical Scheme 

Let us consider the domain: 

1^ = {{x,y) ■■ < X < L,^, yo < y < yo + Ly}. 
On this domain, we address the numerical solution of the following equation: 

i{kxdx + kydy)u + + — = 0, (18) 

where v = i/(x) and /x = /x(x); it is supplemented by the same boundary condition as before on 
{x = Q}: 

icky{kxdy kydx)u + 2kxU = 2kxg^ 

where g is given by Equation ((9|). It is the same problem as in Section [21 except that the coefficients 
u and /X may be functions of x. In the sequel, we consider alternatively the case where /i is a function 
of as a matter of fact, we can take 

/i = /(|n|), where/(t/;) =e-°"''-l, 

with a a positive constant (for a justification of this model, see for example [19] [H] ). 

The interesting problems involve a very small coefficient and it may be necessary to have a 
sufficiently small so that there is no blow-up of the solution. 
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3.1 Description of the Scheme 



Let us set : 



v = I'o + vi with z^o = iiif i^, 

so is a constant and ui a function of x. One discretizes the problem according to a regular grid, 
we denote by 6x, 6y the space step in the two directions and by n and j the indices corresponding 
respectively to x and y; then u{n5x, j5y). 

The numerical method is based on a space marching technique according to the x variable and 
a splitting with respect to this variable. According to Proposition [2l when the value of u" is known, 
it would be possible to evaluate a first intermediate value u'^^^^ by solving on + 6x] the 

following equation: 

{kxdx + kydy)u — i-A_|_u + vqu = 0. 

it would be given by ^{u'''^''') = r{u'')e^-^^'^^^'' (here we denote T = Ty). 

As a matter of fact, in order to have an accurate treatment of the advection term, we prefer to 
perform the following simple splitting : at each space step [x^jx" + one solves succesively 

kxdxU - i^Ax-u + Vi^u = 0, 
kxdxU + kydyU + (i^i + ifi)u = 0. 



3.1.1 Initialization 

For the initial condition, recall that 

g = ie^{kxdy-kydx)u''' + u''', 

where the input data = it|"^Q is a smooth function of the transverse variable y = • x = 
kxy — kyX which values zero around the corner points y = yo and y = yo + Ly, so one can take its 
Fourier transform. 

To determine the boundary value uP of u, we use Formula (fT2]) 



.F(nO) = (19) 



1 + J 1 - 2^ + 2ii^''^'S 

That is to say, {u^)j is obtained by taking the FFT (Fast Fourier Transform) of g, dividing this 

function of ry by the function 1 + ~ + 2iz^*"^ and then taking the IFFT (Inverse Fast 
Fourier Transform) of the result. 

Generally, the input data u*" is a sum of Gaussian functions whose half-height width is in the 
order of a characteristic length Ls which is the typical value of the speckle width (a speckle is a hot 
spot inside the laser beam) and Ls is generally larger than 20 times e. Then one checks that for 
values of e/Lg less than 0.1, the term ieky{kxdy — kydx)u^^ that appears in the previous formula for 
5 is a corrective term and it is possible to take simply g equal to ti*". 
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3.1.2 First stage: Fourier transform 

The first stage is to solve 

kxdxU - i|A_Lti + vqu = 0, (20) 

and we proceed from n" to u"*. Practically, from Proposition O we get immediately : 

^(n'^*) =^(n")e(^-(*'^)+^''^)^^ 
In fact, we have 

ky 2vq 2irje{ri — iv^ky) 



R-i^v)+ivf = , ^ (21) 

" kx{l + ^1 - 2^ + 2ii.of ) kl{l + ^1 - 2^ + 2i^of )2 

Notice that this formula may be used even if is equal to zero, provided that the square root of 
the complex quantity is well defined. 

So, after a FFT on {vJ^), we multiply it by e^'^"*-*''^'''*''^^^'^^ and then apply an inverse FFT. We 
denote ("uj*) the value of the intermediate function, in the cell (n,j). 

3.1.3 Second stage: finite diff'erence scheme 

Boundary conditions on the edges {y = 0} and {y = L} 

It is well known that for this kind of propagation model, the boundary treatment is sensitive; see 
for example [2] for the case of wave equations. In our case the problem is somehow difi'erent since 
there is a privileged direction of propagation: as we use a FFT technique, the key point at each stage 
of the space marching scheme is to force the values of the numerical solution to be negligeable on 
both edges. Therefore we use a damping method which is well known by physicists who address this 
kind of problem [15]. The principle is to introduce in a strip near each edge an artificial absorbing 
coefficient denoted by i?; it decreases progressively on the first five cells near the edge and is very 
large on the edge. More precisely, if i/fj denotes the value of vi in cell {n,j), one replaces i/fj by 

+ Bj where the artificial coefficient Bj is defined by 



Bj = bf3^-^ if j < 5 

= ifJ^ax-i<5 (22) 

= elsewhere, 

with f3 typically in the order of 10 to 100. The numerical tests below (with a characteristic value 
of b in the order of 0.1 to 1) show that this technique leads to get a vanishing value of the solution 
on the edges. One checks on Table [3] that the value of the solution (outside the artificial absorbing 
layers) is almost independant from the choosen values of b and f3. Indeed, near the boundary, the 
main step is the advection one and it is crucial to have a numerical solution which is negligible near 
the boundary cell, in order to avoid a spurious ray to appear on the opposite boundary, due to the 
FFT. Notice that, according to the advection scheme by space marching, the modification in the 
artificial layer at position x"^ has no significant impact on the value outside the artificial layer at 
position x""*"^. 

First order scheme. 
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In this stage, we solve on [x'^,x^ + 5x] the following equation: 

kxdxU + kydyU + i^i(x")u + ifj,u = 0. (23) 

To do this, we use standard finite difference methods. Assume that ky > (the case ky < is 
similar). We consider an upwind method, given that the CFL stability criteria 9 < 1 must be 
checked, where 

Q _ 

kx Sy' 

The initial value is now u"* and we get the final value u^~^^ by setting 

where u^* = + (1 — 9)u^'^ . It is the value of the function on the characteristic line passing 

by (x"~^^,yj); for the first cell, we set = 0. 

For the nonlinear model where the term ^ is replaced by /(|n|), the coefficient /x" has to be 

replaced by f{\u^f\) ■ 

Second order scheme 

When 9 = 1, the previous scheme gives very accurate results, but in real cases it is not possible 
to impose this condition, one has 9 <1 and results are much worse (see Table [2]). We improve the 
numerical scheme when < 1 by using a second order scheme as in all advection problems. To do 
this, we choose a flux-limiter method (see ^7]), with the Van Leer function as limiter (tests prove 
it to be the best one: see Figure [5] and Section [3.3.ip . That is to say, we introduce the function </> 
which depends on the ratio A of the gradient of the function in two neighboring cells: 



We have to solve simultaneously two scalar equations (one for the real and one for the imaginary 
part) with the same flux limiter, so we have to choose one single significant quantity to estimate 



the flux limitor: we choose the energy of the laser, i.e. \u\'^ , and evaluate cf) in terms of |u,p and 



not of \uj\: 

Ul^ |2 

+ \Uj I 

We now replace, in the flrst order scheme, the term derivative in y, uj — uj_^, by Fj — Fj^i where 
the flux Fj is defined as: 

F,=uf + ^il-9)iuf^,-nf)cP{X,). 
The second order scheme is now: 
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3.1.4 Numerical method for two-ray model 

One may also consider a more complex model with two rays crossing each other, with two different 
propagation vectors and (one with positive and one with negative y— component: ky > and 
ky < 0.) To do so, it is necessary to evaluate the nonlinear term /(|n|). Theoretically, the laser 
energy is: 

But we are in the framework of W.K.B. approximation and we do not model the fluctuation of the 
solution at the wavelength level. Hence, the term / has to be taken on a function w corresponding 
to the variation of the index of refraction, which is here the average value of |u| over a wavelength: 



w 



\/|ui|2 + |n2|2_ 



One considers the following model, for p = 1,2: 



ikP • VuP + + iuuP = /(V|ui|2 + |n2|2)uP. 

The first stage of the previous scheme is the same as before : for each ray, we consider Equation 
((201) with its own propagation direction k^ or k^. The interaction between the two rays changes 
only the nonlinear term of the second stage. 

3.2 Properties of the scheme 

3.2.1 Stability 

Let us denote llv""!!!- = Yl P^?/- 

j 

Proposition 4 The numerical first order scheme is monotone decreasing for the t^-norm, i.e. the 
following inequality stands 

VnEN, ||n"||p < ||n"+^||p. (27) 

Moreover, the previous inequality is strict if u ^ 0. 

Proof. 

1. First stage: Fast Fourier Transform 

Let us denote by C the discrete variable associated to rj. On the one hand, since 

,n# ^ /i7i7T(^e(^-(*^)+*'^fe)'^^FFT(u")^ 

and since the FFT conserves the /^-norm, we have: 

||n"#||p = ||e(^-('^)+'^fe)''"FFr(n")||,2. 
On the second hand, the inequality 7^e(ii_(iC)) < implies that 
I e I < 1 



u 
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with an equality iff z/q = 0. We deduce that: 

and conclude: 

||n"#||p<||^."||,2, 
with ||'u"*||p = ||ii"||;2 iff uq = 0. 

2. Second stage: upwind scheme 

For the first order scheme, Relation jMI gives us that: 

n+1 ^ S.^j Sy^^j ^j-^> 2ri,J+^^jj% 

Provided that ^u^* - !^{u]* - u]f^) = ^u^f , we obtain 



^ S^2{^^l±lhJ_^n#, (28) 



Since the modulus of the multiplicative coefficient in the right-hand side is smaller than one, 
this leads to ||u""'"-^||j2 < 1 1 |p • By the triangle inequality: 

\\{u-*)^\y<e\\{uf,)^\\,. + {i-e)\\{uf)^^^^^ 

which concludes the proof. 



In the linear case, that is the case where /i is a data and not a function of the scheme is 
obviously consistent, so Proposition H] implies the convergence of the scheme. 

Concerning the second order scheme modifying the advection step, it is well known (cf |17| ) 
that the effect of this technique with a flux-limiter is to allow small CFL— numbers with a better 
accuracy (than the first order scheme) without generating spurius oscillations. These assertions will 
be confirmed by numerical tests we have performed (see Section [3.3. ip . 

3.2.2 Comparison with the classical Schrodinger equation 

If ky — > 0, Equation (fTSl) reduces to the classical Schrodinger equation, in the case fi = f{\u\) : 

idxU + ^OyyU + iuu- f{\u\)u = 0, (29) 

with a very simple boundary condition (notice that g n*"^) 

ui,=o = n'" (30) 

Proposition 5 If ky ^ 0, the solution given by the numerical scheme converges to the solution of 
the classical Schrodinger problem (flPI )( fg5|) . 



13 



Proof. 

* Initializing. Formula (fT^ used in the scheme shows that 

lim T{u;x = 0) = J^{g), 

SO the boundary condition tends to u\x=o — 9i which is Equation ([30]l . 

* First stage. If ky tends to zero, i.e when the ray tends to be perpendicular to the boundary, 
Formula ([2T]l shows that: 

lim i?_ (iry) + ir?-r^ = —i^ — 'i-^ff' , 

SO m"* given by the first stage is the solution of the classical Schrodinger equation without potential: 

idxU + ^dyyU + iuu = 0, 

which is the limit of the advection-Schrodinger equation. 

* Second stage. It corresponds to a classical discretization of the ordinary differential equation: 

dxU + i^iu + if{\u\)u = 0. 

In other words, the scheme is a classical splitting between dispersion and refraction in the Schrodinger 
equation 

3.3 Numerical results 

Let us recall that the laser energy density is equal to Moreover, the physical meaning of the 
absorption coefficient u is the following: with a constant value of u, if there would be no diffraction 
operator, the laser intensity (integrated on a line orthogonal to the propagation direction) would 
decrease by a factor 1/e^ on a propagation distance equal to 1/u. 

We now give the standard numerical values used for the numerical tests. 

1. For the incoming boundary condition on the edge x = 0, we take a Gaussian of amplitude 1 
centered at a point (0, yo) i-e. u^^ = exp( — (/c^(y — yo) — ^y^)"^ 1^1) with Lg = 2.5 /im; which 
corresponds to the typical half-width of a speckle of a laser beam. 

2. For the incidence angle, we take —45'', then k = (— ■^). 

3. e = 0.05 //m, the wavelength of the laser is 27re « 0.31 ^m. 

4. z/Q = f^i = 5.10"^ ^m~^. Notice that the larger the absorption coefficient, the easier the 
numerical simulation (indeed the laser energy decreases faster with respect to the propagation 
distance). 

5. We take a = 5.10"^. It depends on the electronic density of the plasma: in the vacuum a 
would be null. This size order corresponds either to a dense plasma or to a high laser intensity 
- since we have taken a normalized value of the intensity corresponding to a maximum value 
of u*"^ equal to 1. 

6. For the definition of the boundary layer B, given by ([22]l . we take h = 0.1 and (3 = 50. 



14 




Figure 1: Reference case: 5x = 6y = 0.05, Figure 2: 1st order scheme convergence with 

CFL = 1. Then Lf^c = 59.7, Max(|up) = CFL = 1 as a function of cell size 6x (see 

2.14. Tabled]). 



Number of points 


2« 2' 2« 2^ 2iu 


2^1 


Mesh size 6x = 6y 


1.6 0.8 0.4 0.2 0.1 


0.05 


Error on energy Sj^^ — {u^ ''^\'^\6x6y /\u^'^^\'^ 


46 % 32% 15% 6% 2% 




Focusing distance Lf^c 


82.7 61.4 59.5 59.4 59.9 


59.7 


Error on focusing distance 


38% 2.9% 0.4% 0.6% 0.3% 




Maximum of energy Maxnj{\u^j\''^) 


1.74 2.16 2.13 2.13 2.14 


2.14 


Error on the maximum of energy 


19% 0.7% 0.4% 0.4% 0.07% 





Table 1: Convergence of the scheme, with CFL = 1. The last column represents the fully converged 
reference case u^'^^. 



All our figures represent the laser energy \u\ . 

To be easier to read, our examples are variations with respect to the case defined by the previous 
numerical values of the coefficients and computed with a CFL number 6 equal to 1 (see Figured]). 
With these assumptions, the scheme converges very well as the discretization step decreases (see 
Tabled])- Due to the a coefficient, focusing occurs: the beam focuses and reaches a maximum, then 
decreases. Notice that it may even focus several times for larger values of a. All our comparisons 
are made with this reference case, denoted u^'^\ in the fully converged situation (with mesh size 
6x = 0.05, corresponding to 2^^ points on a domain length Lx = 100.) 

3.3.1 Convergence of the scheme 
Convergence of the first order scheme 

We first take the CFL number equal to 1, which is the case where the first and the second order 
schemes are equivalent. To verify the convergence of the scheme, we have three possible indicators. 
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CFL 


0.5 0.6 0.75 0.875 1 


1 


Error on energy — li^.'^^'^^p n^'^^P 


19 % 17% 14% 9% 2% 


- 


Focusing distance 


43.1 49.1 55.6 48.0 59.9 


59.7 


Error on focusing distance 


28% 18% 7% 19% 0.3% 




Maximum of energy 


1.08 1.18 1.42 1.72 2.14 


2.14 


Error on the maximum of energy 


50% 45% 34% 20% 0.07% 





Table 2: Convergence of the first order scheme with cell size 5y = 0.1 and various CFL. The last 
column represents the fully converged reference case already seen vJ''^^ (with 5y = 0.05). We see 
that the focusing phenomenon is very poorly captured (huge error on the maximum of energy as 
soon as CFL < 1). 

A first indicator is the total energy in the physical domain of interest (that is to say, outside the 
artificial absorbing layer) which is equal to the norm of the energy: we denote it by 

\u\^ = Y,n,j\u^\^5x5y. 

So we compare this quantity to the corresponding one of the fully converged case in the 

two first tables, we give the values of the relative error — '^\'^\5x5y /\u^'^^\'^ for different 

cases. Now, if we want to compare for instance the effects of the variation of the incidence angle, 
two other indicators are more relevant in the framework of the nonlinear model. One is given by the 
focusing distance: we can look for the focusing maximal point i^/oc and we measure the distance 
from Lfoc to the origin of the ray. A last indicator is the maximal value of the energy. These last 
two indicators are quite sensitive. For the nonlinear model, the numerical results are illustrated by 
Figure [2] for the reference case ; the estimates of the indicators are close to the ones of the reference 
case when the spatial step decreases (see Table . 

Thus, we may conclude that when CFL = 1, we reach an accurate result even for 5x = 5y = 0.4, 
and that the focusing phenomenon is very well captured. 

If CFL number decreases, the accuracy becomes bad and even the focusing disappears: see 
Tableland Figures [3] and [H (Of course, if the CFL number is strictly larger than 1, the computed 
solution blows up). 

Convergence of the second order scheme 

We tested three different functions for the fiux limiter: the first one is the Van Leer fiux function 
defined by (f25l) . the second one is a convex combination of Lax-Wendroff and Beam- Warning fiux 
limiter functions, defined by 

r if A < 

(j){X) = { A if < A < 1 (31) 






if 


A < 


A 


if 


< A < 1 


1 


if 


1 < A, 


the Superbee functi 





if 


A < 


2A 


if 


< A < i 


1 


if 


i < A < 1 


A 


if 


1 < A < 2 


2 


if 


2 < A. 



(f){X) = { 1 if i < A < 1 (32) 
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10 20 30 40 50 60 70 80 90 100 

X 

Figure 3: First order scheme with CFL = 0.6, 
5x = 0.1, 5y = 0.17. No focusing observed: 
the convergence of the scheme is poor. 



0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 
CFL 

Figure 4: First order scheme: error on the 
maximum of energy, as a function of CFL (see 
Table [2|). 



We always apply these flux limiter functions at A = \u\'^ and not at the real or imaginary part 
of the solution. As clearly shows Figure [5l it appears that the Van Leer flux function is the one 
which gives the most accurate results. It is particularly clear in terms of the error on the maximum 
of energy : even for small CFL, its estimate is quite accurate contrarily to the first order scheme 
(for CFL = 0.5 , the error is only about 3% with second order scheme but about 50% with first 
order one). 

The smaller the CFL is, the more points are needed to get a correct approximation, as illustrates 
a comparison between Figures [7] and El It is however performed even with 2^ points (that is, with 
6x = 0.2) for CFL = 0.6 for instance, contrarily to the scheme of order one, where no focusing at 
all is observed if CFL = 0.6 even for 5x = 0.1 for instance (see Figure [3]). 

Influence of the artiflcial boundary layer 

In the definition of the artificial absorbing layer B given by ([22]) , we make b and (3 vary, with 
fixed cell sizes Sx = 6y = 0.2 and all the other parameters given by the reference case. We look at 
the value of the total energy for each value of b, f3 (the reference values being b = 0.1, (3 = 50.) The 
results are given in Table [D We check that the sensitivity to the exact values of these coeficients is 
very weak; but it is crucial to have b ^ 0, elseif spurious reflexions may appear on the boundaries. 

3.3.2 Variation of several parameters 
• Variation of the absorption coefficient 

The numerical scheme can also be used with no absorption (z^ = 0), it still works and give good 
results. The repartition of z/q and z/i changes very little the solution, as shows Table [H In each 
case, the reference is taken for i'q = ui = ^. The table shows the results only for the comparison on 
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0.06 




Figure 5: Error on the maximum of energy as 
a function of CFL, for 6y = 0.1, for 3 different 
flux limit ers. 



Figure 6: Incidence of the variation of e on 
the focusing distance (all other parameters as 
in the reference case, except and Ly). 




_0.il ' ' ' ' ' ' ' 1 

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 

8 

X 

Figure 7: CFL = 0.8, second order scheme 
with Van Leer flux limiter: error on the fo- 
cusing phenomenon as a function of the cell 
size Sx. 



_0.il ' ' ' ' ' ' ' 1 

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 

8 

X 

Figure 8: CFL = 0.6, second order scheme 
with Van Leer flux limiter: error on the fo- 
cusing phenomenon as a function of the cell 
size Sx. 
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(3 = 10 /3 = 30 /3 = 50 /? = 100 


b=0 


29% 29% 29% 29% 


b=0.1 


0.08% 0.02% 0.02% 


b=0.2 


0.03% 0.03% 0.05% 0.07% 


b=0.5 


0.08% 0.14% 0.15% 0.16% 


b=l 


0.19% 0.22% 0.23% 0.23% 



Table 3: Incidence of the variation of the boundary layer B on the difference between the total energy 
of each case and the one of the reference case {b = 0.1 and rj = 50): ^j,n\ — \uY^''^\'^\6x6y /\u^'^^\'^ . 
The results of this table show that the influence is negligible, as soon as b is not zero. 





£0 

V 


= 




i^i = 0.3 = 0.5 


^ = 0.7 


£0 _ 

V 


0.9 


£0 

V 


= 1 


reference case: 




















u = 10-3, a = 0.05 


0. 


3% 


0.2% 


0.1% 


0.1% 


0.2 


% 


0. 


3% 


u = 10-3, a = 0.5 


6. 


2% 


5.0% 


2.5% 


2.5% 


5.0 


% 


6. 


2% 


1/ = 10-^ a = 0.05 


0. 


5% 


0.4% 


0.2% 


0.2% 


0.4% 


0. 


5% 


u = 10-^ a = 0.5 


8. 


9% 


7.2% 


3.6% 


3.7% 


7.4% 


9. 


3% 



Table 4: Influence of the repartition between vq and v\ in different cases: percentage of error on 
total energy, defined by — 

the total energy; indeed, the focusing distance remains completely unchanged in any case, and the 
maximum of energy changes by less than 0.3% in the worst case. 

When the absorption coefficient is larger, the problem is easier to solve since the laser energy 
decreases when x increases: for instance in the reference case, if we set v = 10-^ instead of = IQ-^, 
the ray is rapidly totally absorbed, and no focusing is observed. 

The influence of the repartition between and v\ increases with a, as shows Table [H 

• Variation of the incidence angle 

To test whether the scheme is accurate for various angles, we make it vary from 5° to 70°, all the 
other parameters being constant: see Table [5l We check that the indicators for the focusing distance 
and the maximum of energy are well estimated, since they depend very few on the incidence angle. 

• Variation of e 

If all other coefficients are fixed, the larger e becomes, the more important the diffusion phe- 
nomenon is (and the larger the domain must be to obtain a converging solution), and, in the 
nonlinear case, the smaller the focusing distance becomes. A limit value of e is experienced, above 
which no focusing phenomenon (for the nonlinear equation) is observed. In our reference case for in- 
stance, the limit is around e = 0.17, see Figure [6l but this limit depends of course on all parameters, 
especially a and v. 

From a physical point of view, all our asymptotic analysis is built on the assumption e = o(l) : 
else, our equation is no more a valid approximation of the envelope of Helmholtz equation, given 
by ([I]). Hence, we have to assume e << 1 : larger values are meaningless. 
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Incidence angle 


5" 


30° 


45" 


60" 


70" 




0.23 


0.16 


0.2 


0.16 


0.02 




0.02 


0.1 


0.2 


0.27 


0.06 


Maximum of energy 


2.17 


2.16 


2.13 


2.10 


1.99 


Error on the maximum of energy 


1.5% 


0.8% 


0.43% 


2.2% 


9.7% 


Focusing distance 


59.2 


59.7 


59.35 


59.9 


60.2 


Error on the Focusing distance 


0.9% 


0.01% 


0.6% 


0.34% 


0.96% 



Table 5: Variation of the incidence angle: influence on the focusing distance and on the maximum 
of energy. As usual, the errors refer to the fully- converged reference case. 



• Variation of a. 

The parameter a represents a nonlinear effect, and induces autofocusing and filamentation of 
the beam. The larger it is, the more accurate the focusing phenomenon becomes, as illustrated in 
Figure M 

It could be interesting to evaluate the value of a for which a focusing phenomenon appears: in 
our reference case, it is for a > 0.02. On the other hand, one may check that if a is large enough, 
several focusing points appear and a breaking of the beam occurs (see Figure [TTI). This phenomenon 
depends of course also on the absorption coefficient u and on the diffusion coefficient e. 



3.3.3 Remark on artificial damping 

We wish to check now that there is no artificial damping due to the numerical scheme; in other 
words, that in the second stage the decrease of the P— norm of the solution has the right value. 
Using the notations of Section [321 this right value is given by the equality: 

Going back to Equation (f28ll . we can write it under the form (assuming no artificial boundary layer: 
B, = 0) 

- 1 + a + ib ' 

where we set a = -^^i j and b = j^-l^^j- Since the characteristic value of the coefficient a is 10~^ 
(or smaller) and, in the worst case, the characteristic value of /i is in the order of 1, so that we can 
choose ^ to have b small, we see that 

\- = 1 - 4a- -r + o(a^), 

which is very close to the right value e~^"' = 1 — 4a + o(a^). The only damping may then come 

from the fact that ^ l^g^P may be significantly smaller than Yl due to a large difference 

j ^ j 
between nj* and To check this numerically, we test the case u = : Figure [TU] shows that 

even in a difficult case with a large a = 1.5, the global energy Hu^Hp is conserved. 
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0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 
a 

Figure 9: Influence of a on the maximum 
of energy (obtained in the focusing phe- 
nomenon). Standard hypothesis. The autofo- 
cusing, which is a nonlinear effect, is more 
significant when a increases. 



-0.012 1 ^ ^ ^ ^ ^ 1 

10 20 30 40 50 60 

Figure 10: a = 1.5, = : we define the 
energy i?" = This picture shows 

(i?" — E^)/E'^ as a function of = n5x : the 
energy E'^ decreases by less than 2% during 
the whole trajectory. 



3.3.4 Two-ray model 

We have also performed computations for the two-ray model which is described above at Section 
13.1.41 using two functions and u^; an illustration is given by Figure [121 The interaction between 
the rays is only given by the nonlinear term f{w) with w"^ = + as above. To analyse its 
exact influence, one can compare the result given by the previous model with the two-ray interaction 
and the result given by a simple superposition of two independant rays (obtained with the one-ray 
model). One may see then that the energy becomes larger with the two-ray interaction: on the case 
of Figure [T^ for instance, Max{\u^\'^ -\- |n^P) = 12.3 instead of 10.6 if the rays do not interact. 

4 Extension to a Time-Dependent Interaction Model 

We now address a model where a tilted paraxial equation is coupled with a hydrodynamic model in 
order to study filamentation. Under the hypothesis of a small incidence angle, this model has been 
extensively used by physicists for a long time and it is also addressed in |3],[3],[l0] for example and 
the references therein (for a derivation of this model, see [TH] for example). 

4.1 The Model and the Numerical Method 

Modeling of the plasma. 

By taking the critical density (depending only on the laser wave length) as a reference density, 
one defines a non-dimension electron density N = N{t, x) ; so the plasma may be characterized 
only by this quantity, the plasma velocity U = U(t,x) and the electron density Te(t,x). 
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10 20 30 40 50 60 



Figure 11: a = 1.5, z/ = : high focusing. 
One observes a breaking of the beam in three 
sub-beams. 






□ 20 40 BO ao 100 133 140 160 

Figure 12: 2 beams crossing with incidence 
angles ±30^^, a = 0.05, and L = 5 for the 
initial gaussian functions. 



Then, the simplest model is the following one. The pressure P = P{N,Te) is assumed to be a 
smooth function of the density and of the electron temperature Tg (which is assumed to be a 
very smooth fixed function of the position x ), for example P{N,Te) may be the sum of two terms 
equal to A^^ and NTe up to multiplicative constants. Then one considers the following barotropic 
Euler system: 

^iV + V(iVU) = 0, (33) 

— (A^U) +V(A^UU) +V(P(A^,re)) = -iV7pV|^'|2. (34) 

The term 7pV|^'p corresponds to a ponderomotive force due to a laser pressure (the coefficient 
7p is a constant depending only on the ion species). 
Modeling of the laser beam. 

The laser field = ^'(t, x) is a solution to the following frequency wave equation (which is of 
Schrodinger type): 

Id 1 

2i-—'if + —A'^ + ko(l-N)^ + iiy'''^ = 0, (35) 
c at kq 

where the real coefficient is related to the absorption of the laser intensity by the plasma and c 
the light speed. 

Assume that the mean value of the plasma density is quite constant and denoted N^, so we set: 

N{^)=Nm + 6N{^), 

where 6N is small with respect to 1. Then one can make the paraxial approximation ; that is to say 
the laser beam is now characterized by the space and time envelope of the electric field U = h({t, x) 
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and we set: 



^'(t, x) = U{t, x)e'*^"^-'', where K = y^l - iV^k. 



Therefore, if one sets e = i by the same procedure as mentioned in the introduction, one 

checks that U satisfies: 

h5N^^ IdU 



Vl - Nmi^.VU + -^.iU) + i—U - 4— + = 0. (36) 

2 2 2 c ai 



It is necessary to supplement equation (|36]) with the same boundary condition as in the model 
of section 1 (and with an initial condition). 
Numerical method. 

We consider a mesh of finite difference type as above. The numerical treatment of the barotropic 
Euler system ((33]) ([Ml) is a classical one, we have chosen a Lagrange- Euler method, see [3] for details. 
To deal with (f36]) . according to the large value of the speed of light, one must perform a time inplicit 
discretization. So at each time step, one solves firstly the Euler system with a ponderomotive force 
evaluated with the previous value of Secondly, using the obtained values of N and of 5N, one 
has to solve ([36]) ; if u*™ and u denote the values of the field U at the beginning and the end of 
time step, one searches u solution to: 

ik.Vu + iuu + ^{Aiu)- fiu= ' (37) 

2 CVl — Nrn ot 

where we have set: 

ko5N 11 1 o 

u = — , u = — — H V . 

2^/1 -Nj cVl - Nm St 2V1 - iVm 

That is exactly the equation studied in section 3, but a right hand side term has been added. So 
the numerical method is the same as described above ; the only modification is the adding of the 
right hand side term in the transport stage. Notice that the index of refraction (1 — N) is equal to 
(l-2e^)(l-iV^). 

>From a practical point of view, the numerical method for ([36]) has been implemented in a 
parallel way in the HERA plateform for plasma hydrodynamics in 2D and in 3D; the parallel solver 
and the domain decomposition techniques are the same as the ones detailed in [3]. 



4.2 Numerical Results 

Recall that from a practical point of view, in the transverse profile of a laser beam, one distinguishes 
a lot of small hot spots, called speckles, whose intensity is very large compared to the mean intensity 
of the beam. The shape of each individual speckle is a Gaussian function whose width is about 
a few micrometers. We present here the results of a 2D numerical simulation. One addresses a 
simulation box which is 600 /im long and 300 fim wide, the laser propagates with an incidence angle 
of 19*^. The incoming boundary condition a = a{y) is independent of time and mimics a laser beam 
whose width is equal to 40/im with five speckles ; each speckle is modeled by a centered Gaussian 
function h and is characterized by a random phase Cfc, that is to say a{y) = 'E^^^akh{y — yfc)e*^'=, 
where the are random and the are close to each other. The plasma has an initial density 
equal to Nm = 0.15 and the temperature is equal to 35. 10^ Kelvin. The mesh consists of 4 millions 
of cells and the time step is in the order of 0.1 picosecond (it is determined at each time step by 
the Courant-Friedrichs-Levy condition related to the sound speed of the plasma). The initial value 
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Figure 13: Snapshot of the laser intensity at the time 2.6 ps, 3.9 ps, 5.3 ps and 6.6 ps ( from the top-left to 
the bottom- right). 

of the laser intensity is zero, the plasma is progressively grabed by the ponderomotive force and on 
Figure [iHl we have plotted the laser intensity at different times. At the first snapshot (at time 2.6 
ps), the plasma is not grabed enough, so the value of /i is small; the autofocusing effect is very low 
but not negligible: instead of five different speckles at the incoming boundary one notices only four 
speckels at the rear side (one of the four has a larger intensity) and a little spreading of the beam 
may be observed. At the second snapshot, the position of the four speckles has changed and the 
plasma is more grabed - since the energy density is larger in one speckle. On the two last snapshots, 
we may check that the spreading of the beam at the rear side of the simulation box becomes larger 
when the time increases. Moreover the configuration is not stationary, this situation is characteristic 
of the so-called filamentation instability. 

Conclusion 

A mathematical analysis has lead to an analytical form of the solution of the tilted paraxial equation 
in the simple case where the refraction index and the absorption coefficients are constant. After- 
wards, we proposed a numerical method for solving the initial problem which uses the previous 
analytical form. The scheme has the property to yield a classical scheme when incidence angle 
becomes zero and the equation reduces to the classical paraxial one. The numerical method is illus- 
trated by some results on toy problems. We have also given extensions of this model, which have 
enlarged the capability of our plateform HERA for laser propagation in a plasma (see |lj and [HI 
for examples of simulations performed with HERA). This numerical method may be also extended 
in the case where the unit vector K depends slowly on the one-dimension spatial variable x.n, for 
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instance if one has to deal with an equation of the following type 

iK.Vu + i-(V.K)u + - fiu + ivu = 0, on V. 

2 2ko 

The paraxial equation in a tilted frame may be also considered in a first region where the plasma 
density is slowly varying with respect to the spatial variable and coupled with another model in a 
neighbor region where the plasma density is strongly varying: in that region the laser is no more 
characterized by the time-space envelope of the fast oscillating electric field but by the wave equation 
(f35ll (see [6], for results obtained in HERA with this model). For simulating such a physical tilted 
beam, a classical paraxial model without accounting for the incidence angle would lead to search a 
the solution which would be highly oscillating with respect to the space variable and therefore to 
increase dramatically the mesh size to get accurate results. 
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